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The critical point of lattice QCD on the /i-T plane 
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t^j- | Abstract: We propose a method to study lattice QCD at finite temperature and chem- 

ical potential. We compare it with direct results and with the Glasgow method by using 
n/=4 QCD at Im(/i)7^0. We locate the critical endpoint (E) of QCD on the Re(/x)- 
T plane. In this study we use nf— 2+1 dynamical staggered quarks with semi-realistic 
masses on L t =A lattices. 



QCD at finite T and [i is of fundamental importance, since it describes relevant features 
of particle physics in the early universe, in neutron stars and in heavy ion collisions. Exten- 
sive experimental work has been done with heavy ion collisions at CERN and Brookhaven 
to explore the fi-T phase diagram. Note, that past, present and future heavy ion exper- 
iments with always higher and higher energies produce states closer and closer to the T 
axis of the \i-T diagram. It is a long-standing non-perturbative question, whether a critical 
point exists on the [i-T plane, and particularly how to tell its location theoretically (T|. 

Let us discuss first the /i=0 case. Universal arguments pi and lattice simulations 
Q indicate that in a hypothetical QCD with a strange (s) quark mass (m s ) as small as 
the up (u) and down (d) quark masses (m u d) there would be a first order finite T phase 
transition. The other extreme case (nj=2) with light u/d quarks but with an infinitely 
large m s there would be no phase transition only an analytical crossover. Note, that 
observables change rapidly during a crossover, but no singularities appear. Between the 
two extremes there is a critical strange mass {m c s ) at which one has a second order finite T 
phase transition. Staggered lattice results on Lj=4 lattices with two light quarks and m s 
around the transition T (n/=2+l) indicated E| that m c s is about half of the physical m s . 
Thus, in the real world we probably have a crossover. 
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Returning to a non- vanishing /i, one realizes that arguments based on a variety of 
models (see e.g. || || |J) predict a first order finite T phase transition at large /i. Combining 
the fi = and large fi informations an interesting picture emerges on the \i-T plane. For 
the physical m s the first order phase transitions at large fi should be connected with the 
crossover on the fi = axis. This suggests that the phase diagram features a critical 
endpoint E (with chemical potential (j>e and temperature Te), at which the line of first 
order phase transitions (fi > fiE and T < Te) ends pi]. At this point the phase transition 
is of second order and long wavelength fluctuations appear, which results in (see e.g. @) 
consequences, similar to critical opalescence. Passing close enough to (/j,e,Te) one expects 
simultaneous appearance of signatures which exhibit nonmonotonic dependence on the 
control parameters [||, since one can miss the critical point on either of two sides. 

The location of E is an unambiguous, non-perturbative prediction of QCD. No ab 
initio, lattice analysis based on QCD was done to locate the endpoint. Crude models with 
m s = oo were used (e.g. suggesting that fiE ~ 700 MeV, which should be smaller for 
finite m s . The goal of our exploratory work is to show how to locate the endpoint by a 
lattice QCD calculation. We use full QCD with dynamical n/=2+l staggered quarks. 

QCD at finite /i can be formulated on the lattice ||; however, standard Monte-Carlo 
techniques can not be used. The reason is that for Re(^)^0 the determinant of the Eu- 
clidean Dirac operator is complex. This fact spoils any importance sampling method. 

An attractive approach to alleviate the problem is the "Glasgow method" (see e.g. 
Ref. in which the partition function (Z) is expanded in powers of exp(/x/T) by using 
an ensemble of configurations weighted by the n=0 action. After collecting more than 20 
million configurations only unphysical results were obtained. The reason is that the /i=0 
ensemble does not overlap enough with the finite density states of interest . 
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We propose a method to reduce the 
overlap problem and determine the phase 
diagram in the /i-T plane (for details see 
jl2]l ). The idea is to produce an ensemble 
of QCD configurations at /U=0 and at T c . 
Then we determine the Boltzmann weights 
1 14 ] of these configurations at fj, ^ and at 
T lowered to the transition temperatures 
at this non-vanishing [i. Since transition 
configurations are reweighted to transition 
configurations a much better overlap can 
be observed than by reweighting pure had- 
ronic configurations to transition ones [|10|]. 
After illustrating the applicability of the 

method we locate the critical point of QCD. (Multi-dimensional reweighting was successful 
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Figure 1: ipip as a function of lm(\x), for di- 
rect results ( squares ), our technique ( crosses ) and 
Glasgow-type reweighting (dots). 



for determining the endpoint of the hot electroweak plasma [13] e.g. on 4D lattices.) 

Let us study a generic system of fermions ip and bosons 4>, where the fermion Lagrange 
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density is ipM{(p)ip. Integrating over the Grassmann fields we get: 

Z(a) = J V<Pexp[-S bos (a, <f>)] det M(<f>, a), (1) 

where a denotes a set of parameters of the Lagrangian. In the case of staggered QCD a 
consists of (3, m q and [/,. For some choice of the parameters a=ao importance sampling 
can be done (e.g. for Re(/x)=0). Rewriting the above equation we get 

Z(a) = j V(f) exp[- S bos (a ,(j))] det M(4>, a ) 

f r ^ / ,s „ , ,m det M(4>, a) 1 , . 

|exp[-^(a, 0) + S bos (a , 0)] — ^ — } j . (2) 

We treat the curly bracket as an observable (measured on each configuration) and the rest 
as the measure. Changing only one parameter of the ensemble generated at «o provides an 
accurate value for some observables only for high statistics. This is ensured by rare fluctu- 
ations as the mismatched measure occasionally sampled the regions where the integrand is 
large. This is the overlap problem. Having several parameters the set ao can be adjusted 
to get a better overlap than obtained by varying only one parameter. 

The basic idea of the method as applied to dynamical QCD can be summarized as 
follows. We study the system at Re(//)=0 around its transition point. Using a Glasgow- 
type technique we calculate the determinants for each configuration for a set of [i, which, 
similarly to the Ferrenberg-Swendsen method 14] , can be used for reweighting. The average 



plaquette values can be used to perform an additional reweighting in (3. Since transition 
configurations were reweighted to transition ones a much better overlap can be observed 
than by reweighting pure hadronic configurations to transition ones as done by the Glasgow- 



type techniques (moving along the transition line was also suggested by Ref. |17|| ). 

We have directly tested these ideas in n/=4 QCD with m g =0.05 dynamical staggered 
quarks. We first collected 1200 independent V=4-6 3 configurations at Re(^)=Im(//)=0 and 
some values and used the Glasgow-reweighting and also our technique to study Re(/i)=0, 
Im(^)/0. At Re(/u)=0, lm(/i)^0 direct simulations are possible. After performing these 
direct simulations as well, a clear comparison can be done. Figure 1 shows the predictions 
of the three methods for the average quark condensates at /3=5.085 as a function of Im(/x). 
The predictions of our method agree with the direct results, whereas for larger Im(^) the 
predictions of the Glasgow method are by several standard deviations off. Based on these 
experiences we expect that our method can be applied at Re(/i)^0. 

In QCD with nj staggered quarks one should change the determinants to their nj/A 
power in our two equations. Importance sampling works also in this case at some (3 and at 
Re(^)=0. Since detM is complex an additional problem arises, one should choose among 
the possible Riemann-sheets of the fractional power in eq. (||). This can be done by using 
the fact that at fj, = [a w the ratio of the determinants is 1 and the ratio should be a 
continuous function of [i. However, the continuity can only be ensured if the analytical 
dependence of the determinant on fj, is known |l^] (the idea goes back to a method of pq]). 
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Figure 2: 
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ImfPy ) as a function of the chemical 



In the following we keep \i real and look 
for the zeros of Z for complex (3. At a 
first order phase transition the free energy 
oc log Z{(5) is non-analytic. A phase transi- 
tion appears only in the V^ co limit, but 
not in a finite V. Nevertheless, Z has zeros 
at finite V, generating the non-analyticity 
of the free energy, the Lee- Yang zeros . 
These are at complex values of the parame- 
ters, in our case at complex j3. For a system 
with a first order transition these zeros ap- 
proach the real axis in the oo limit (de- 
tailed analysis suggests \jV scaling). This 

co limit generates the non-analyticity of the free energy. For a system with crossover 
Z is analytic, and the zeros do not approach the real axis in the V^ co limit. 

At T/0 we used L t =4, L s =4,6,8 lat- 
tices. T=0 runs were done on 10 3 - 16 lat- 
tices. m Uj d=0.025 and m s =0.2 were our 
bare quark masses. At T ^ we deter- 
mined the complex valued Lee- Yang zeros, 
(3q, for different V-s as a function of \x. 
Their V— > co limit was given by a (3q(V) = 
/?o° + C/V extrapolation. We used 14000, 
3600 and 840 configurations on L s =4,6 and 
8 lattices, respectively. Im(/3^°) is shown 
on Figure 2 as a function of ll. For small 
fj,-s the extrapolated Im(/3^°) is inconsis- 
tent with a vanishing value, and predicts a 
crossover. Increasing ll the value of Im(/3(J°) 
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Figure 3: 
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The T-li diagram. Direct results 
are given with errorbars. Dotted line shows the 
crossover, solid line the first order transition. The 
box gives the uncertainties of the endpoint. 



decreases, thus the transition becomes con- 
sistent with a first order phase transition. (Note, that systematic overshooting is a finite 
V effect.) The statistical error was determined by jackknife samples of the total L s = 4, 6 
and 8 ensembles. Our primary result is n e nd = 

0.375(20). 

To set the physical scale we used a weighted average of Ro, m p and \fa. Note, that 
(including systematics due to finite V) we have (Rq -m n ) = 0.73(6), which is at least twice, 
m U: d is at least four times as large as the physical values. 

Figure 3 shows the phase diagram in physical units, thus T as a function of lib, the 
baryonic chemical potential (which is three times larger then the quark chemical potential). 
The endpoint is at T E = 160 ± 3.5 MeV, ll e = 725 ± 35 MeV. At Lt B =0 we obtained 
T c = 172 ± 3 MeV. 

We proposed a method -an overlap improving multi-parameter reweighting technique- 
to numerically study non-zero fi and determine the phase diagram in the T-fj, plane. Our 
method is applicable to any number of Wilson or staggered quarks. As a direct test 
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we showed that for lm(/i)^0 the predictions of our method are in complete agreement 
with the direct simulations, whereas the Glasgow method suffers from the well-known 
overlap problem. We studied the fi-T phase diagram of QCD with dynamical nf =2+1 
quarks. Using our method we obtained Te=160±3.5 MeV and //£=725±35 MeV for the 
endpoint. Though (j,e is too large to be studied at RHIC or LHC, the endpoint would 
probably move closer to the /i=0 axis when the quark masses get reduced. At ^=0 we 
obtained T c =172±3 MeV. Clearly, more work is needed to get the final values. One has to 
extrapolate in the R-algorithm and to the thermodynamic, chiral and continuum limits. 

This work was partially supported by Hung. Sci. grants No. OTKA-T34980/T29803/- 
T22929/M28413/OM-MU-708/IKTA111/NIIF. This work was in part based on the MILC 
collaboration's public lattice gauge theory code: |http://physics.indiana.edu/~sg/milc.htmi 
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